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A  conference  in  numerical  linear  algebra  was  held  at  the  Mathematical 
Sciences  Research  Institute  at  the  University  of  California  at  Berkeley  on 
October  17.  1993.  The  occasion  was  the  00th  birthdays  of  Professors  Beres- 
ford  Parlett  and  Yelvel  Kalian.  There  were  approximately  100  attendees.  In 
addition  to  Denmtel.  Prof.  James  Bunch  of  UC'  San  Diego  and  Dr.  Horst 
Simon  of  NASA-Ames  were  co-organizers.  Here  is  the  list  of  speakers: 

•  James  Bunch,  University  of  California.  San  Diego.  "Three  Decades  of 
Numerical  Linear  Algebra  at  Berkeley". 

•  G.  W.  Stewart,  University  of  Maryland,  "On  the  Perturbation  of  Ma¬ 
trix  Factorizations". 

•  John  Reid.  Rutherford  Appleton  Laboratory.  England.  ’'Taking  Ad¬ 
vantage  of  Sparsity  within  2x2  Pivots  When  Solving  Symmetric  Indef¬ 
inite  Sets  of  Linear  Equations". 

•  Horst  Simon.  NASA  Ames.  "Spectral  Algorithms- A  New  Approach  to 
Some  Discrete  Optimization  Problems  in  Scientific  Computing". 

•  Larry  Nazareth.  Washington  State  University,  "The  Newton  and  Cauchy 
Perspectives  on  Computational  Nonlinear  Optimization". 
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•  James  Ortega.  University  of  Virginia.  "Solution  of  Nonlinear  Poisson- 
type  Equations". 

•  Anne  Greenbaum,  Courant  Institute,  NYU,  "Matrices  that  Generate 
the  Same  Krylov  Residua!  Spaces". 

•  Bahrain  Nour-Omid,  Scopus,  Berkeley.  "Ordered  Modified  Gram-Schmidt". 

•  David  Scott.  Intel.  "A  High  Performance  Out-of  Core  Dense  Equation 
Solver  for  the  Intel  Parallel  Supercomputer". 

•  Scott  Baden.  University  of  California.  San  Diego.  The  Role  of  Heuris¬ 
tics  in  Parallel  Computation  of  Scientific  Problems". 

•  Peter  Tang.  Argonne  National  Laboratory.  "Recent  Advances  in  Rank- 
Revealing  QR  Factorization". 

•  James  Dernmel.  University  of  California.  Berkeley.  "Recent  progress 
in  parallel  algorithms  for  the  nonsymmetric  eigenproblern". 

The  speakers  anc'  others  submitted  papers  for  a  special  issue  of  the  Jour¬ 
nal  of  Numerical  Linear  Algebra  with  Applications  (JNLAA)  dedicated  to 
Parlett  and  Kalian.  These  papers  were  sent  out  for  the  standard  refereeing 
process.  Most  have  been  accepted,  and  several  are  still  out  for  review  and 
revision.  We  attach  a  list  of  these  papers,  abstracts  and  their  status  below. 
When  the  refereeing  process  is  complete  and  the  journal  published,  we  will 
forward  copies  to  ONR  as  requested. 


Titles . txt 


JNLAA 


lCIAL  issue  dedicated  to  paelett  and  kahan 


AUT 

HORS 

(1) 

Demme  1 , 

Higharo,  and  Schreiber 

(2) 

Eisner , 

He,  and  Mehrmann 

(3)  Fierro  and  Bunch 

(4)  Ikramov 

(5)  Z.  Jia 

(6)  Paige,  Parlett,  and 

van  der  Vorst 

(7)  Van  Huff  el  ar.d  Park 
(3)  Arbenz  and  Golub 

(3)  Borges,  Frezza,  and  Gragg 

(10)  Druskin  and  Knizhnerman 

(11)  C-T  Pan  and  Peter  Tang 

(12)  Nour-Omid,  Dunbar  and  Woodbury 

(13)  Eatterson 

(14)  P.-C  Li 

(15)  Z.  Bai 

(16)  Nazareth 

(17)  Edelman  and  Mascarenhas 
(13)  Dan  S  zy 1 d 

(19)  Walden,  Karl son,  and  Sun 

(2  0;  o.  Stewart 


Block  LU  Factorization 

Minimization  of  the  norm,  the  norm  of  the 
inverse  and  the  condition  number  of  a  matrix 
by  completion 

Orthogonal  Projection  and  Total  Least  Squares 
On  the  Computational  Aftermath  of  Matrix 


Consimilarity  and 

Con  comm: 

utativity 

Relations  II 

The  Convergence  of 

Krylov 

Subspace 

Methods 

Approximate  Soluti 

or.s  and 

Eigenval 

ue  Bounds  frc 

Krylov  Subspaces 

Efficient  Algorithms  for  Bordered  Band  Matrices 

Matrix  Shapes  Invariant  under  the  QR  Algorithm 

Some  Inverse  Problems  for  Jacobi  and  Arrow  Matr 

Krylov  Subspace  Approximation  of  Eigenpairs  anc 
Matrix  Functions  in  Exact  and  Computer  Arithmet 

Bounds  on  Singular  Values  Revealed  by  QR 
Factorizations 

Ordered  Modified  Gram- Schmidt.  Or  thogonali  zatior 

Dynamical  Analysis  of  Numerical  Systems 

Solving  Secular  Equations  Stably  and  Efficient! 

Progress  in  the  Numerical  Solution  of  the 
Nonsymmetric  Eigenvalue  Problem 

Trust  Regions  based  on  Conic  Functions  in 
Linear  and  Nonlinear  Programming 

On  Parlett 's  Matrix  Norm  Inequality 

Regions  of  Convergence  of  the  Rayleigh 
Quotient  Iteration  Method 

Optimal  Backward  Perturbation  Bounds  the 

I ’ "ear  Least  Squares  Problem 

On  the  Solution  of  Block  Hessenberg  Systems 


Apr  i 


.59) 


3 


BLOCK  LU  FACTORIZATION 


JAMES  W.  DEMMEL 

Computer  Science  Division  and  Mathematics  Department 
University  of  California 
Berkeley 
CA  94720 
U.S.A. 

{na.deueltna-net  .oral. gov). 


NICHOLAS  J.  HIGHAM 
Nuffield  Science  Research  Fellow 
Department  of  Mathematics 
University  of  Manchester 
Manchester 
MIS  9PL 
UK 

(na.nhighaa®na-net .  oral . gov). 


ROBERT  S.  SCHREIBER 

RIACS 

NASA  Ames  Research  Center 
Moffett  Field 
CA  9 4  OS  5 
U.S.A. 

(na .  schre  ibex  *na-net .  oral .  go*). 


ABSTRACT 

Many  of  the  currently  popular  “block  algorithms"  are  scalar  algorithms  in 
which  the  operations  have  been  grouped  and  reordered  into  matrix  operations. 
One  genuine  block  algorithm  in  practical  use  is  block  LU  factorization,  and 
this  has  recently  been  shown  by  Demmel  and  Higham  to  be  unstable  in  general. 
It  is  shown  here  that  block  LU  factorization  is  stable  if  A  is  block  diagonally 
dominant  by  columns.  Moreover,  for  a  general  matrix  the  level  of  instability 
in  block  LU  factorization  can  be  bounded  in  terms  of  the  condition  number 
k(A)  and  the  growth  factor  for  Gaussian  elimination  without  pivoting.  A 
consequence  is  that  block  LU  factorization  is  stable  for  a  matrix  A  that  is 
symmetric  positive  definite  or  point  diagonally  dominant  by  rows  or  columns 
as  long  as  A  is  well-conditioned. 

Keywords-,  block  algorithm,  LAPACK,  level  3  BLAS,  iterative  refinement,  LU 
factorization,  backward  error  analysis,  block  diagonal  dominance. 

AMS(MOS)  subject  classifications,  primary  65FG5,  65F25,  65G05. 


1  Introduction 

Block  methods  in  matrix  computations  are  widely  recognised  as  being  able  to 
achieve  high  performance  on  modern  vector  and  parallel  computers.  Their  per¬ 
formance  benefits  have  been  investigated  by  various  authors  over  the  last  decade 
(see,  for  example,  [11,  14,  15J),  and  in  particular  by  the  developers  of  LAPACK  [1], 
The  rise  to  prominence  of  block  methods  has  been  accompanied  by  the  development 
of  the  level  3  Basic  Linear  Algebra  Subprograms  (BLAS3) — a  set  of  specifications  of 
Fortran  primitives  for  various  types  of  matrix  multiplication,  together  with  solution 
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Minimization  of  the  norm,  the  norm  of 
the  inverse  and  the  condition  number  of 
a  matrix  by  completion 

L.  Eisner  C.  He  V.  Mehrmann 

Fakultat  fur  Mathematik 
Universitat  Bielefeld 
Postfach  8640 
D-4800  Bielefeld,  FRG 


Abstract 

We  study  the  problem  of  minimizing  the  norm,  the  norm  of  the  inverse  and 
the  condition  number  with  respect  to  the  spectral  norm,  when  a  submatrix 
of  a  matrix  can  be  chosen  arbitrarily.  For  the  norm  minimization  problem  we 
give  a  different  proof  than  that  given  by  Davis/Kahan/Weinberger.  This  new 
approach  cam  then  also  be  used  to  characterize  the  completions  that  minimize 
the  norm  of  the  inverse.  For  the  problem  of  optimizing  the  condition  number 
we  give  a  partial  result. 

Keywords:  condition  number,  norm  of  a  matrix,  matrix  completion,  dilation 
theory,  robust  regularization  of  descriptor  systems 


1  Introduction 


We  study  the  following  optimization  problem:  Given  integers  n,m,N  >  n,m  and 
matrices  A  €  Cn-m,B  €  Cn<N~m,  C  €  CN~n'm,  find  X  €  Cs~nS~m  such  that  the 


matrix 


(1.1) 


satisfies 

cond(W(X))  =  minireC„_.,„_..{|lW(Z)||j|W(Z)-Mi} 

(1.2) 

=  minZ6C*_».,,_f.{cond(W(Z))}. 
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Orthogonal  Projection  and  Total  Least  Squares 


(3) 


Ricardo  D.  Fierro  and  James  R.  Bunch 

Department  of  Mathematics 
University  of  California,  San  Diego 

La  Jolla,  CA  92093 

Dedicated  to  Beresford  Parlett  and  William  Kahan 
on  the  occasion  of  their  60th  birthdays. 


ABSTRACT.  When  the  overdetermined  system  of  linear  equations  AX  as  B  has  no 
solution,  compatibility  may  be  restored  by  an  orthogonal  projection  method.  The 
idea  is  to  determine  an  orthogonal  projection  matrix  P  by  some  method  M  such 
that  [A  B ]  =  P[A  B],  and  AX  as  B  is  compatible.  Since  both  A  and  B  are  allowed 
to  be  perturbed,  denote  by  Xm  the  minimum  norm  total  least  squares  solution  by 
method  M  (the  TLS-M  solution)  to  AX  =  B.  In  this  paper  conditions  for  com¬ 
patibility  of  the  lower  rank  approximation  and  subspace  properties  of  A  in  relation 
to  the  nearest  rank-fc  matrix  to  A  are  discussed.  We  find  upper  and  lower  bounds 
for  the  difference  in  the  TLS-M  solution  Xm  and  the  SVD-based  TLS-SVD  solu¬ 
tion  Xsvd  and  also  provide  a  perturbation  result  for  the  SVD-based  TLS  method. 
These  results  suggest  a  new  algorithm  for  computing  a  total  least  squares  solution 
based  on  a  rank  revealing  QR  factorization  and  subspace  refinement.  Numerical 
simulations  are  included  to  illustrate  the  conclusions. 
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APPROXIMATE  SOLUTIONS  AND  EIGENVALUE  BOUNDS 
FROM  KRYLOV  SUBSPACES 


CHRIS  C.  PAIGE 
Computer  Science 
McGill  University 
Montreal,  HSA  SA7 
Canada 


BERESFORD  N.  PARLETT 
Mathematics  Department 
University  of  California 
Berkeley,  CA  94720 


HENK  A.  VAN  DER  VORST 
Mathematical  Institute 
University  of  Utrecht 
Budapestlaan  6 
Utrecht,  the  Netherlands 


The  first  and  third  authors  dedicate  this  to  Bcresford  Parlett  and  Velvet 
Kahan  —  not  only  in  recognition  of  their  exceptional  abilities  in 
general,  and  contributions  to  this  topic  in  particular,  but  also  for  the 
friendship,  ideas  and  encouragement  they  have  given  us,  and  above  all 
for  the  impeccable  character  and  style  they  have  both  consistently 

exhibited. 


ABSTRACT 

Approximations  to  the  solution  of  a  large  sparse  symmetric  system  of  equa¬ 
tions  are  considered.  The  conjugate  gradient  and  minimum  residual  approxi¬ 
mations  are  studied  without  reference  to  their  computation.  Several  different 
bases  for  the  associated  Krylov  subspace  are  used,  including  the  usual  Lanc- 
zos  basis.  The  zeros  of  the  iteration  polynomial  for  the  minimum  residual 
approximation  are  given  the  name  harmonic  Ritz  values  and  are  character¬ 
ized  in  several  ways  and,  in  addition,  attractive  convergence  properties  are 
established.  The  connection  of  these  harmonic  Ritz  values  to  Lehmann’s  op¬ 
timal  intervals  for  eigenvalues  of  the  original  matrix  appears  to  be  new. 

Keywords:  Krylov  subspace;  Lanczos  process;  symmetric  matrix;  conjugate 
gradients;  minimum  residual;  Lehmann  intervals. 
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EFFICIENT  REDUCTION  ALGORITHMS 
FOR  BORDERED  BAND  MATRICES 
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Dedicated  to  Beresford  Parlett  and  William  Kahan  on  the  occasion  of  their  60th  birthdays. 


Abstract 

Various  plane  rotation  patterns  are  presented,  which  provide  stable  algorithms  for 
reducing  a  6-band  matrix  bordered  by  p  rows  and/or  columns  to  (6  +  p)-band  form. 
These  schemes  generalize  previously  presented  0(N 5)  reduction  algorithms  for  ma¬ 
trices  of  order  N,  6  =  1.  and  p  =  1  to  the  reduction  of  more  general  6-band, 
p-bordered  matrices  where  6  >  1  and  p  >  1.  Moreover,  by  splitting  the  matrix  into 
two  similarly  structured  submatrices  and  chasing  nonzeros  to  the  corners  in  two  di¬ 
rections,  the  newly  proposed  patterns  reduce  the  number  of  required  rotations  and 
hence  the  computational  cost  by  one  half  compared  to  the  other  existing  one-way 
chasing  algorithms.  Symmetric,  as  well  as  more  general  matrices,  are  considered. 
An  example  of  the  first  type  is  the  symmetric  arrowhead  matrix  that  arises  in  solv¬ 
ing  inverse  eigenvalue  problems.  Examples  of  the  second  type  are  found  in  updating 
the  singular  value  decomposition  (SVD)  and  the  partial  SVD. 


Keywords:  arrowhead  matrix,  band  matrix,  inverse  eigenvalue  problem,  Givens 
rotations,  singular  value  decomposition,  updating. 


1  Introduction 

A  matrix  of  order  N  is  said  to  be  b-band  or  b-diagonal  if  it  consists  of  6  diagonals. 
If  p  trailing  rows  are  attached,  then  the  resulting  matrix  A(.v+p)x^r  is  called  a 
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Some  Inverse  Eigenproblems  for  Jacobi  and  Arrow  Matrices 


CARLOS  F.  BORGES 
Code  Ma/Bc 

Naval  Postgraduate  School 
Monterey,  CA  9S94S 

W.B.  GRAGG 


RUGGERO  FREZZA 
DEI,  Univ.  di  Padova 
via  Gradenigo  6/ A 
35131  PADOVA  ■  ITALY 


Code  Ma/Gr 

Naval  Postgraduate  School 
Monterey,  CA  93943 


ABSTRACT 

We  consider  ihe  problem  of  reconstructing  Jacobi  matrices  and  real  symmet¬ 
ric  arrow  matrices  from  two  eigenpairs.  Algorithms  for  solving  these  inverse 
problems  are  presented.  We  show  that  there  are  reasonable  conditions  un¬ 
der  which  this  reconstruction  is  always  possible.  Moreover,  it  is  seen  that  in 
certain  cases  reconstruction  can  proceed  with  little  or  no  cancellation.  The 
algorithm  is  particularly  elegant  for  the  tridiagonal  matrix  associated  with  a 
bidiagonal  singular  value  decomposition. 

Keywords:  Jacobi  matrix,  Arrow  matrix,  inverse  problem. 


1  Introduction 

We  consider  the  problem  of  reconstructing  Jacobi  matrices  and  real  symmetric  arrow 
matrices  from  two  eigenpairs.  The  algorithms  we  present  for  solving  these  inverse 
problems  are  simple,  and  useful  for  constructing  test  matrices  for  eigenproblems. 
The  algorithm  for  reconstructing  Jacobi  matrices  was  applied  to  the  problem  of 
model  identification  of  reciprocal  stochastic  processes  in  [3]. 


2  Jacobi  matrices 

Let  T  be  an  unreduced  real  symmetric  tridiagonal  matrix  (i.e.  a  Jacobi  matrix) 


A 

A  ih 


(23) 


with  A  >  0  for  i  =  1,2, ...,  n  -  1.  We  use  the  notation  introduced  in  [13]  and  let 
UST(n)  denote  the  set  of  n  x  n  real  unreduced  symmetric  tridiagonal  matrices, 
and  let  UST+(n)  denote  that  subset  of  IJST(n)  with  positive  A- 
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abstract 

The  modified  Gram-Schmidt  algorithm  is  used  to  orthogonalize  a  vector,  r.  against 
a  set  of  orthogonal  vectors,  u  and  v.  By  means  of  an  error  analysis,  it  is  shown  that  the 
order  in  which  the  vectors  u  and  v  are  processed  can  have  a  significant  effect  on  the  state 
of  orthogonality  of  the  final  set  of  vectors.  Alt  ordered  modified  Gram-Schmidt  algoritlun 
is  presented  in  wliich  the  order  of  orthogonalization  is  determined  by  computing  the 
magnitude  of  the  components  of  r  along  the  vectors  of  an  orthogonal  basis  and  the 
orthogonalization  proceeds  in  the  order  of  descending  component  magnitude.  The  benefit 
of  the  ordered  algoritlun  is  an  improved  orthogonality  stale.  Examples  involving  the  use 
of  the  Arnold!  algoritlun,  in  which  the  modified  Gram-Sdunidl  algoritlun  is  used  to 
obtain  a  set  of  orthogonal  basis  vectors,  demonstrate  the  effect  of  the  order  in  which 
the  ortliogonalizations  are  performed  and  that  faster  convergence  can  result  when  the 
ordered  algoritlun  is  used  to  obtain  an  orthogonal  set  of  basis  vectors. 

Keywords:  Gram-Sclunidt  orthogonalization,  unsynunetric  matrices,  Krylov  subspace, 
Arnold!  algorithm 


1.  Introduction 

In  this  paper  the  problem  of  orthogonalizing  a  given  vector,  r,  against  a  vector 
basis,  {qj ,q;, . . .  ,qm},  which  is  orthogonal  to  working  precision,  is  addressed.  In 
practical  situations,  the  intent  is  to  augment  the  basis  with  an  orthogonalized  and 
normalized  r.  Tne  modified  Gram-Schmidt  (MGS)  algorithm  is  typically  used  for 
this  purpose.  The  focus  of  the  paper  is  on  the  order  in  which  the  vectors  q,  are 
processed  in  MGS.  Although  in  exact  arithmetic  the  result  is  independent  of  order 
in  which  ortliogonalizations  are  performed,  the  following  illustrative  example  and 
error  analysis  show  that  in  finite  precision  the  order  of  orthogonalizations  can  have 
a  significant  effect  on  the  state  of  orthogonality  of  the  final  set  of  vectors. 
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DYNAMICAL  ANALYSIS  OF  NUMERICAL  SYSTEMS 
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STEVE  BATTERSOS 

Department  of  Mathematics  and  Computer  Science 
Emory  University 
Atlanta,  GA  <0322  USA 
sb@mathes.em0ry.edu 


ABSTRACT 

For  many  years  techniques  from  numerical  analysis  have  been  applied  fruitfully  ;o  the  study 
of  dynamical  systems.  In  this  paper  it  is  shown  that  the  theory  of  dynamical  systems  may  be  applied 
to  certain  computational  problems  In  particular  the  question  of  global  convergence  of  vanous  QR 
algorithms  can  be  reduced  to  the  study  of  certain  vector  iterations  derived  from  Sdiur  forms  of 
matrices.  The  technique  is  illustrated  in  determining  the  convergence  behavior  of  normal  Hessenberg 
matrices  under  the  Francis  and  multishift  QR  iterations. 


The  purpose  of  this  paper  is  lo  demonstrate  how  techniques  from  dynamical  sys¬ 
tems  can  be  applied  to  certain  problems  in  numerical  analysis.  Iterative  algoriihms  such 
as  Newton  s  method,  the  power  method,  and  the  QR  algorithms  may  be  viewed  as 
dynamical  systems.  In  this  setting  a  variety  of  dynamical  teenniques  address  numerical 
issues  of  convergence.  We  will  place  previous  Rayleigh  quotient  iteration  work  of  ParJett 
and  Kahan  into  this  framework  and  explain  how  it  generalizes  to  our  results  on  the  global 
convergence  of  Francis  shifted  QR. 

The  author  was  originally  trained  in  the  theory  of  dynamical  systems  and  later 
developed  an  interest  in  eigenvalue  computation.  This  paper  is  directed  at  a  numerical 
analysis  audience.  It  is  the  author's  intention  to  persuade  the  audience  of  a  fruitful  inter¬ 
disciplinary  connection  between  his  interests. 

1.  The  Iterative  Algorithm-Dynamical  Systems  Paradigm 
This  section  answers  the  following  three  questions: 

(i)  What  is  a  dynamical  system? 

(ii)  Describe  a  central  problem  in  dynamical  systems. 

(iii)  What  computational  questions  can  be  placed  in  the  setting  of  question  ii? 

While  there  are  several  possibilities  we  define  a  (discrete)  dynamical  system  to  be 
a  map  from  a  topological  space  M  into  itself.  In  most  cases  the  space  will  have  additional 
structure  (e.g.  manifold)  and  the  map  have  some  degree  of  smoothness.  For  the  present 
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THE  CONVERGENCE  OF  KRYLOV 
SUBSPACE  METHODS  FOR  LARGE 
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Abstract 

The  convergence  of  both  FOM  and  GMRES  for  solving  large  unsymmetric 
linear  systems  is  still  a  hard  problem.  Saad  and  some  others  have  partially 
solved  it,  mainly  for  the  case  where  the  coefficient  matrix  A  is  diagonalizable 
and  its  spectrum  lies  in  the  open  right  (left)  half  plane.  In  this  paper,  we 
will  focus  on  the  convergence  problem  of  both  FOM  and  GMRES  in  the  case 
where  the  coefficient  matrix  A  is  defective.  When  the  spectrum  of  A  either  lies 
in  the  open  right  (left)  half  plane,  namely,  (A  +  AH)f 2  is  positive  (negative) 
definite,  or  is  on  the  real  axis,  we  establish  the  related  theoretical  error  bounds 
and  reveal  the  intrinsic  relationships  between  the  convergence  speed  and  the 
spectrum  of  A.  Besides,  we  make  the  two  more  notes  for  the  case  where 
A  is  diagonalizable.  The  results  show  tha!  FOM  and  GMRES  ale  likely  1 
converge  slowly  when  either  A  is  defective  oi  the  Jardan  basis  ol  A  is  ill 
conditioned.  Meanwhile,  we  point  out  that  most  results  derived  in  this  paper 
are  also  suitable  for  the  convergence  analysis  of  some  other  methc  s,  eg.  GCR 
etc. 


Keywords:  unsymmetric  linear  systems,  Krylov  subspace,  the  Chebyshev 
polynomials,  diagonalizable,  defective,  derivatives. 
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Krylov  subspace  approximation  of  eigenpairs  and  matrix 

FUNCTIONS  IN  EXACT  AND  COMPUTER  ARITHMETIC 

C  SO  C 

Vladimir  DRUSKIN  ,  Leonid  KNIZHNERMAN 

C  *o 

Schlumfcerger  Doll  Research,  Old  Quarry  Road,  Ridgefield, 
CT  06377-4108,  the  U.  S.  A.  Telephone  C  203}  431-54-71,  telefax 
C  203D  438-38-19,  E-mail  DRUSKI  N€>sdr  .  si  b.  com. 

C  »«o 

Central  Geophysical  Expedition,  Russia,  123298,  Moscow, 
Narodncgo  Opolcheniya  St.  ,  house  40,  build.  3.  Telephone 
C  0955  192-81-56,  telex  411758  ATLAS  SU,  E-mai 1 

KNI ZHNEROsei smos. msk . su. 


Abstract . 

Many  researchers  are  now  working  on  computing  the 

product  of  a  ma+*~ix  function  and  a  vector,  using 

approximations  in  the  Krylov  subspace.  We  review  our  results 
on  the  analysis  of  one  interpretation  of  that  approach  for 
symmetric  matrices,  which  we  call  Spectral  Lanczos 

Decomposition  Method  C SLDMD  . 

There  has  been  proved  a  general  convergence  estimate, 
relating  SLDM  error  bounds  with  those  obtained  through 
approximation  of  the  matrix  function  by  a  piece  of  its 
Chebyshev  series.  Thus,  we  arrived  at  effective  estimates 
for  matrix  functions,  arising  when  solving  parabolic, 
hyperbolic  and  elliptic  partial  derivative  equations.  We 
concentrate  on  the  parabolic  case,  where  we  obtain 

estimates,  indicative  of  the  super convergence  of  SLDM.  For 
this  case  a  symbiosis  of  SLDM  and  splitting  is  also 
considered  and  some  numerical  results  are  presented. 
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Further,  we  implement  our  general  estimates  for  getting 
convergence  bounds  of  Lanczos  approximations  to  eigenvalues 
in  the  internal  part  of  the  spectrum  that,  unlike 
Kaniel-Saad  estimates,  are  independent  of  the  set  of 
eigenvalues,  located  between  the  required  one  and  the 
nearest  spectrum  hound. 

We  consider  an  extension  of  our  general  estimate  for  the 
case  of  the  simple  Lanczos  method  C without 
reorthogonal i zationD  in  the  computer  arithmetic,  which  shows 
that  for  a  moderate  dimension  of  the  Krylov  subspace  the 
results,  proved  for  the  exact  arithmetic,  are  stable  to 
roundoi f . 

We  also  touch  upon  the  Arnoldi  method. 

§  1.  Introduction 

This  paper  is  a  review  of  the  authors’  results  on 
analysis  of  Krylov  approximations  to  a  matrix  function  by  a 
vector  and  of  the  Lanczos  method,  partially  published  in  the 
former  Soviet  Union. 

Let  A  be  a  symmetric  nxn  matrix.  It  is  the  Lanczos 
method  that  has  become  accepted  as  a  powerful  tool  for 
finding  eigenpairs  of  A,  at  least,  if  A  is  large  and  sparse. 
However,  the  Lanczos  method  has  got  another  application. 
Namely,  let  /  be  a  function,  defined  on  the  spectral 
interval  of  A,  and  £>eRn.  Let  us  consider  computation  of  the 
vector 

\i=/CA><p.  Cl  5 

Solving  systems  of  linear  equations  is  a  problem  of  this 
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Solving  Secular  Equations 
Stably  and  Efficiently 

Ren-Cang  Li 

Department  of  Mathematics 
University  of  California  at  Berkeley 
Berkeley,  California  94720 

October  16,  1992 


Dedicated  to  B.  N.  Parlett  and  H7.  Kahan  on  the  occasion  of  their  60th  birthdays 

Abstract 

The  way  described  in  [2]  to  solve  secular  equations,  which  is  the 
kernel  of  the  divide  and  conquer  method  for  solving  symmetric  tridi¬ 
agonal  matrix  eigensystem,  is  analyzed.  New  and  more  efficient  ap¬ 
proaches  are  proposed.  Some  practical  issues  concerning  implementa¬ 
tions  of  these  approaches  are  tackled  in  detail. 
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Progress  in  the  Numerical  Solution  of 
the  Xonsymmetric  Eigenvalue  Problem* 

Dcdicaltd  to  IP.  Kahan  and  B.  N.  Parlttt  on  iht  occasion  of  their  60lh  birthdays 

Zhaojun  Bai  t 
February  22.  1993 


Abstract 

With  the  growing  demands  from  disciplinary  and  interdisciplinary  fields  of  science 
and  engineering  for  the  numerical  solution  of  the  nonsymmetric  eigenvalue  problem, 
competitive  new  techniques  have  been  developed  for  solving  the  problem.  In  this  paper, 
we  examine  the  start-of-the-art  of  the  algorithmic  techniques  and  the  software  scene  for 
the  problem.  Some  current  developments  are  also  outlined. 


1  Introduction 

Over  several  years  working  on  the  LAPACK  project  [2],  and  on  algorithm  and  software 
development  of  the  nonsymmetric  eigenvalue  problem  and  communication  with  a  variety  of 
users  who  work  in  diverse  fields  involving  scientific  computing,  the  author  has  seen  a  growing 
demand  for  the  numerical  solution  of  the  nonsymmetric  eigenvalue  problems.  Meanwhile, 
in  numerical  analysis  community,  since  Parlett's  exploratory  review  paper  entitled  “The 
Software  Scene  in  the  Extraction  of  Eigenvalues  from  Sparse  Matrices”  [34]  nearly  one 
decade  ago,  and  with  the  successful  development  of  the  symmetric  eigenvalue  problem, 
many  new  numerical  methods  and  analysis  have  been  developed  for  the  nonsymmetric 
eigenproblem.  The  aim  of  this  essay  is  to  review  the  origins  of  the  problem^nd  the  progress 
of  the  numerical  techniques  over  the  past  decade,  and  to  share  our  view  ana  expertise  within 
scientific  computing  community. 

The  survey  is  by  no  means  complete.  One  reason  for  this  is  that  relevant  articles  may¬ 
be  found  scattered  throughout  the  scientific  and  engineering  literature,  and  the  task  of 
tracking  them  all  down  is  impossibly  large.  The  author  apologizes  for  the  ignorance  of 
some  important  contributions  to  the  problem  that  are  not  mentioned  here.  A  new  book 
by  Saad  [42]  is  an  elegant  source  for  studying  the  start-of-the-art  in  large  eigenproblem 
techniques.  This  review  will  only  focus  on  the  nonsymmetric  eigenvalue  problem  in  the 
aspects  of  its  origins,  algorithmic  techniques,  software  scene  and  work  in  progress. 

As  defined  by  Parlett  [34]  one  decade  ago,  there  are  two  different  user  groups  for  the 
eigenproblem.  One  is  called  the  intensive  user  group  and  the  other  called  the  sporadic 

‘This  work  was  supported  in  part  by  NSF  grant  ASC-9102963  and  in  part  by  the  Applied  and  Compula- 
tiona)  Mathematics  Program,  Defense  Advanced  Research  Projects  Agency,  under  contract  DM28E04120. 

1  Department  of  Mathematics,  University  of  Kentucky,  Lexington,  KY  40506,  U.S.A,  e-mai]  address: 
baiCas.uky.edu 
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Trust  Regions  Based  on  Conic  Functions  in 
Linear  and  Nonlinear  Programming 

J.L.  Nazareth* 

October.  1992 

Dedicated  to  Professors  B.N.  Parlett  and  W.  Kahan 
on  the  occasion  of  their  6CTth  birthdays 


Abstract 

An  optimization  method  is  developed  based  on  ellipsoidal  trust 
regions  that  are  defined  by  conic  functions.  It  provides  a  powerful 
unifying  theory  from  which  can  be  derived  a  variety  of  interesting  and 
potentially  useful  optimization  algorithms,  in  particular,  conjugate- 
gradient-like  algorithms  for  nonlinear  minimization  and  Karmarkar- 
like  interior-point  algorithms  for  linear  programming. 


1  Introduction 

We  propose  a  new  optimization  method  based  on  the  use  of  trust  regions 
(More  [7],  Ye  [13].  Gonzaga  [5])  defined  by  conic  functions  (Davidon  [1]). 
Our  trust  regions  are  ellipsoidal  with  a  distinguished  axis  of  orientation. 
We  investigate  their  use  in  the  setting  of  nonlinear  minimization  where  they 
lead  to  conjugate-gradient-like  algorithms,  and  in  the  setting  of  linear  pro¬ 
gramming  where  they  lead  to  Karmarkar-like  interior-point  algorithms. 

Conic- based  ellipsoidal  trust  regions,  a  natural  generalization  of  quadratic- 
based  trust  regions,  are  thus  shown  to  provide  both  a  powerful  unifying 
theory  and  a  means  for  formulating  a  variety  of  interesting  and  potentially 
useful  optimization  algorithms. 

‘Professor,  Department  of  PiiTe  and  Applied  Mathematics,  Washington  State  Univer¬ 
sity,  Pullman,  WA  99164-3113  and  Affiliate  Professor.  Department  of  Applied  Mathemat¬ 
ics,  University  of  Washington.  Seattle,  WA  98195.  email:  nazarethQwsumath.bitnet  and 
nazarethQamath.washington.edu 
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OPTIMAL  BACKWARD  PERTURBATION 

BOUNDS 

FOR  THE  LINEAR  LEAST  SQUARES  PROBLEM 

Bertil  Walden,  Rune  Karlson,  and  Ji-guang  Sun  * 

Department  of  Mathematics,  Linkoping  University,  S-58I  88  Linkoping,  Sweden 


Abstract 

Let  A  be  an  m  x  n  matrix,  b  be  an  m-vector,  and  i  be  a  purported  solution  to 
the  problem  of  minimizing  jjfc  -  Ax j];.  We  consider  the  following  open  problem: 
to  find  the  smallest  perturbation  E  of  A  such  that  the  vector  x  exactly  minimizes 
|)6  -  (A  +  E)x IJj.  This  problem  is  completely  solved  when  E  is  measured  in  the 
Frobenius  norm.  When  instead  using  the  spectral  norm  of  E,  easily  computed 
upper  and  lower  bounds  are  given,  and  the  optimum  is  found  under  certain 
conditions. 

AMS  Subject  Classifications:  65F20 

Key  words:  linear  least  squares,  backward  perturbations 


Throughout  this  paper  we  will  use  the  following  notation.  C”*xn  denotes  the 
set  of  complex  m  x  n  matrices,  Cn  =  C"xl.  I  is  the  identity  matrix,  and  0  is  the 
null  mattix.  AT  and  AH  stand  for  the  transpose  and  conjugate  transpose  of  a  ma¬ 
trix  A ,  respectively.  A *  denotes  the  Moore- Penrose  inverse  of  A.  j|  H2  denotes  the 
Euclidean  vector  norm  and  the  spectra]  norm,  and  |j  j(p  denotes  the  Frobenius  norm. 

Let  A  6  Cmxn,6  €  Cm,  and  let  i  be  a  purported  solution  to  the  problem  of  min¬ 
imizing  ||6  -  Ax H2-  Stewart  [3]  discovered  two  perturbations  £  of  A  such  that  the 
vector  x  exactly  minimizes  ||6-  {A  4  £)x||2.  Recently,  Sun  [5]  found  a  perturbation 
E.  (see  below  (10))  and  showed  that  !|£.(|f  is  an  improved  backward  perturbation 
bound,  but  it  has  not  been  proved  that  the  scalar  ||£.|[f  is  the  optimal  backward 
perturbation  bound  in  the  meaning  of  the  Frobenius  norm. 

This  paper  gives  a  solution  to  the  following  open  problem:  to  find  the  smallest 
perturbation  £  of  A  such  that  the  vector  z  exactly  minimizes  j(6  -  (A  +  £)z||2  (ref. 
[1],  [3],  (4,  p.160-163]).  We  use  €  to  denote  the  set  of  all  perturbation  matrices  £ 
of  A  such  that  x  is  an  exact  solution  to  minx  jjh  -  (A  4  £)*!ta-  The  expression  of 

‘The  research  of  this  author  was  supported  by  the  Swedish  Natural  Science  Research  Council. 
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On  the  Solution  of 
Block  Hessenberg  Systems 

G.  W.  Stewart 

Dedicated  to  Velvel  Kahan  and  Beresford  Parlett 
ABSTRACT 

This  paper  describes  a  divide-and-conquer  strategy  for  solving  block 
Hessenberg  systems.  For  dense  matrices  the  method  is  as  efficient  as 
Gaussian  elimination;  however,  because  it  works  almost  entirely  with 
the  original  blocks,  it  is  much  more  efficient  for  sparse  matrices  or  ma¬ 
trices  whose  blocks  can  be  generated  on  the  fly.  For  Toeplitz  matrices, 
the  algorithm  can  be  combined  with  the  fast  Fourier  transform  to  give 
a  new  superfast  algorithm. 


1.  Introduction 

This  paper  was  motivated  by  the  attempt  to  find  the  steady-state  of  Markov 
chains  of  types  M/G/l  and  G/M/l  (see  [8]  for  definitions  and  further  details). 
The  matrix  of  transition  probabilities  of  a  chain  of  type  M/G/l  has  the  block 
upper  Hessenberg  form 

^  Bq  B\  B2  B3  B4 

Co  Ai  A2  A3  A4 

0  Ao  A\  A2  A3  • 

M=  0  0  Ao  Ax  A2  ■ 

0  0  0  Aq  Ai 

The  transition  matrices  of  chains  of  type  G/M/l  are  block  lower  Hessenberg 
matrices,  formed  in  analogy  with  (1.1).  Note  that  after  the  first  row  and  column 
of  (1.1)  are  deleted,  the  matrix  becomes  block  Toeplitz. 

We  shall  return  to  these  chains  at  the  end  of  this  paper.  For  now  we  are 
going  to  consider  the  more  general  problem  of  solving  the  block  upper  Hessenberg 


(1.1) 
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Bounds  on  Singular  Values  Revealed  by  QR  Factorizations 

C'.-T.  Pan 

Department  of  Mathematical  Sciences 
Northern  Illinois  University 
De  Kalb,  IL  60115-2SSS  “ 

Ping  Tak  Peter  Tang* 

Mathematics  and  Computer  Science  Division 
Argonne  National  Laboratory 
9700  South  Cass  Ave. 

Argonne,  IL  GO-J 39-4S0 1 


Abstract 


We  introduce  a  pair  of  dual  concepts:  pivoted  blocks  and  reverse 
pivoted  blocks.  These  blocks  are  the  outcome  of  a  special  column 
pivoting  strategy  in  QR  factorization.  Our  main  result  is  that  under 
such  a  column  pivoting  strategy,  the  QR  factorization  of  a  given 
matrix  can  give  tight  estimates  on  any  two  a  priori-chosen  consecutive 
singular  values  of  that  matrix.  In  particular,  a  rank-revealing  QR 
factorization  is  guaranteed  when  the  two  chosen  ronsecutive  singular 
values  straddle  a  gap  in  the  singular  value  spectrum  that  gives  rise 
to  the  rank  degeneracy  of  the  given  matrix.  The  pivoting  strategy, 
called  cyclic  pivoting,  can  be  viewed  as  a  generalization  of  Golub’s 
column  pivoting  and  Stewart’s  reverse  column  pivoting.  Numerical 
experiments  confirm  the  tight  estimates  that  our  theory  asserts. 

AMS  classification:  65F30,  15A23,  15A42,  15A15. 

Key  words  and  phrases:  Singular  value  decomposition,  rank- 
revealing  QR  factorization,  cyclic  column  pivoting. 


‘This  author  was  supported  by  the  Applied  Mathematical  Sciences  subprogram  of  the  Office  of  Energy 
Research,  U.S.  Department  of  Energy,  under  Contract  W-31-109-Eng-38. 
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On  Parlett’s  matrix  norm  inequality 


A 


Alan  Edelman* 
Walter  F.  Mascarenhas  * 

February,  1993 


Abstract 

We  show  that  a  certain  matrix  norm  ratio  studied  by  Parlett  has  a  supremum  that 
is  0(\/n)  when  the  chosen  norm  is  the  Frobenius  norm,  while  it  is  O(logn)  for  the  two 
norm.  This  ratio  arises  in  Parlett’s  analysis  of  the  Cholesky  decomposition.  We  are 
happy  to  dedicate  this  note  to  our  friends  Beresford  and  Velvel  on  the  occasion  of  their 
sixtieth  birthdays. 


1  Introduction 


This  note  bounds  the  quantity 


x(«) 


=  sup  t(U)  — 


m 

\\UT +  U +UTU\\' 


defined  by  Parlett  [3]  as  the  supremum  over  all  (non-zero)  upper  triangular  matrices  U  with 
diagonal  entries  u,-,-  >  —  1.  In  particular,  we  study  xf{ «)  and  Xz(n)>  where  the  norms  in  \(n) 
are  taken  to  be  the  Frobenius  norm  and  the  two  norm  respectively.  The  quantity  x(n)  arises 
in  Parlett’s  analysis  of  the  Cholesky  decomposition.  The  term  UTU  in  the  denominator 
would  be  neglected  by  first  order  perturbation  theory  but,  according  to  Parlett,  it  actually 
helps  in  the  analysis. 

An  equivalent  expression  for  x(n)  is 


*(n)  =  “p/w  =  jrrjnh’ 

where  the  supremum  is  over  the  set  of  n  by  n  upper  triangular  matrices  N  with  non-negative 
diagonals,  excluding  the  identity  matrix. 

‘Department  of  Mathematics  and  Lawrence  Berkeley  Laboratory,  University  of  California,  Berkeley,  CA 
94720,  edelman@math.berkeley.edu.  Supported  by  the  Applied  Mathematical  Sciences  subprogram  of  the 
Office  of  Energy  Research,  U.S.  Department  of  Energy  under  Contract  DE-AC03-76SF00098. 

’Dept,  de  Matematica,  Universidade  Estadual  de Campinas,  Caixa  Postal  6065,  Barao  Geraldo,  Campinas 
SP  13081,  Brazil,  walterm@dcc.unicamp.br 
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REGIONS  OF  CONVERGENCE  OF  THE  RAYLEIGH  QUOTIENT 

ITERATION  METHOD 

RICARDO  D.  PANTAZIS'  and  DANIEL  B.  SZYLD1 


Dedicated  to  Professoi's  IF.  Kahan  and  B.  N.  Parlett  on  the  occasion  of  their 

sixtieth  birthdays. 

Abstract.  Parlett  and  Kalian  have  shown  in  19CS  that  for  almost  any  initial  vector  in  the  unit 
sphere,  the  Rayleigh  quotient  iteration  method  converges  to  some  eigenvalue.  In  this  paper,  the 
regions  of  the  unit  sphere  which  include  all  possible  initial  vectors  converging  to  a  specific  eigenvalue 
are  studied.  It  is  shown  that  when  the  matrix  is  shifted,  or  multipied  by  a  scalar,  the  regions 
do  not  change.  In  the  case  of  three  dimensions,  these  regions  are  completely  characterized.  It  is 
shown  experimentally  that,  for  the  three-dimensional  case,  the  probability  that  an  initial  vector  will 
lead  to  convergence  to  an  interior  eigenvalue  is  larger  than  the  probability  that  it  converges  to  an 
extreme  one.  It  is  conjectured  that  the  same  is  true  for  matrices  of  any  order.  Experiments  in  higher 
dimensions  are  presented  which  conform  with  the  conjecture. 

Key  words.  Eigenvalues,  eigenvectors,  symmetric  matrices,  Rayleigh  quotient  iteration,  con¬ 
vergence,  basins  of  attraction. 

AMS(MOS)  subject  classification.  C5F15 


1.  Introduction.  We  study  global  convergence  properties  of  the  Rayleigh  quo¬ 
tient  iteration  method  for  the  solution  of  the  symmetric  generalized  eigenvalue  prob¬ 
lem 

(1)  Ax  —  XBx. 

'  Department  of  Computer  Science,  Duke  University,  Durham,  North  Carolina  27706-2591 ,  USA 
(rdp'Scs.duke.ed  u). 
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(szy)d®euclid. math. temple.edu).  The  work  of  this  author  was  supported  by  a  Summer  Research 
Fellowship  at  Temple  University  and  by  the  National  Science  Foundation  grant  DMS-9201728. 
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